#ifndef DYN_PLANARBILLIARD_H
#define DYN_PLANARBILLIARD_H

#include <armadillo>

#include "dyn_planarBilliard_arcLengthFunction.h"
#include "dyn_planarBilliard_nextImpactFunction.h"
#include "dyn_planarBilliard_nextImpactDiff.h"
#include "../utilities/dyn_functor.h"
#include "../utilities/dyn_multifunctor.h"
#include "../ode/dyn_ode.h"
#include "../quad/dyn_quad.h"
#include "../roots/dyn_roots_hunting.h"
#include "../roots/dyn_roots_newtonRaphson.h"
#include "../roots/dyn_roots_multiDim_sta.h"

using namespace arma;

/*! \author Joey Dumont <joey.dumont@gmail.com>
 *
 * \date 2012-10-03
 *
 * \brief Provides functionality to draw trajectories in a 2D billiard.
 *
 * We provide functionality to compute the trajectories of the
 * 2D billiard given the boundary curve as a function of the angle,
 * i.e. \f$r=r(\varphi)\f$.
 *
 * Given a set of initial conditions, the trajectories are propagated
 * in space and specularly reflected at the boundary. We keep track
 * of the trajectory by using the \f$\varphi\f$, the position on the
 * curve, and \f$\alpha\f$, the angular between the tangent to the
 * curve at that point and the trajectory vector. It is also possible to
 * obtain the (normalized) arc length,
 *  \f[
 *      S(\varphi) = \frac{1}{L}\int_0^\varphi\sqrt{r^2(\varphi')+(d/d\varphi')^2}d\varphi'
 *  \f]
 * by numerical quadrature.  The user supplies the quadrature algorithm. The momentum
 * \f$p\equiv\cos\alpha\f$ is
 * also computed. This can be used to draw phase space diagrams.
 *
 * The user can also access the jacobian matrices \f$\mathbf{M}_{10}\f$, which are
 * computed automatically. The radius of curvature, given by
 *  \f[
 *      \rho(\varphi)=\frac{(r^2+r'^2)^\frac{3}{2}}{r^2+2r'^2-rr''},
 *  \f]
 * is computed as an aside the jacobian matrices.
 */

class PlanarBilliard
{
public:
   /*! Constructor needs the boundary curve and its first and second derivatives. */
    PlanarBilliard(Functor& _boundary,
                   Functor& _diffBoundary,
                   Functor& _diff2Boundary,
                   Quad& _quadrature,
                   colvec _initCond,
                   double _tolQuad,
                   double _tolBis,
                   double _tolNew,
                   int _numberBounces);

    /*! @name Accessor Functions */
    //@{
    Functor& getBoundary(){return boundary;}
    Functor& getDiffBoundary(){return diffBoundary;}
    Functor& getDiff2Boundary(){return diff2Boundary;}
    Quad& getQuadrature(){return quadrature;}
    colvec getInitialConditions(){return initCond;}
    double getToleranceQuadrature(){return tolQuad;}
    double getToleranceBisection(){return tolBis;}
    double getToleranceNewton(){return tolNew;}
    int getNumberBounces(){return numberBounces;}

    void setBoundary(Functor& _boundary){boundary=_boundary;}
    void setDiffBoundary(Functor& _diffBoundary){diffBoundary=_diffBoundary;}
    void setDiff2Boundary(Functor& _diff2Boundary){diff2Boundary=_diff2Boundary;}
    //void setQuadrature(Quad& _quadrature){quadrature=_quadrature;}
    void setInitialConditions(colvec _initCond){initCond=_initCond;}
    void setToleranceQuadrature(double _tolQuad){tolQuad=_tolQuad;}
    void setToleranceBisection(double _tolBis){tolBis=_tolBis;}
    void setToleranceNewton(double _tolNew){tolNew=_tolNew;}
    void setNumberBounces(int _numberBounces){numberBounces=_numberBounces;}
    //@}

    /*! Member functions. */
    //@{
    /*! Computes the arc length from theta=0.0 to theta=phi. */
    double computeArcLength(double phi);
    /*! Computes the momentum for given alpha. */
    double computeMomentum(double alpha);
    /*! Gives the next impact with the boundary */
    colvec computeNextImpact(colvec previousImpact);
    /*! Computes the whole trajectory. */
    mat computeTrajectory();
    /*! Finds periodic orbits of given period with given accuraccy. */
    mat huntUPOs(double lambdaInit, double tolSTA, int period, int maxIterations);
    /*! Computes the radius of curvature at a certain angle. */
    double computeRadiusOfCurvature(double phi);
    /*! Computes the jacobian matrix at a particular point. */
    mat computeJacobianMatrix(colvec impact);
    //@}

    /*! We define our own modulus function. */
    double user_fmod(double x, double y)
    {
        return x - floor(x/y)*y;
    }

    /*! Compute STA correction for billiard. */
    colvec computeSTAcorrection(double lambda, mat reflection, colvec previousStep, int period);

protected:
    /*! @name Member variables. */
    //@{
    Functor& boundary;
    Functor& diffBoundary;
    Functor& diff2Boundary;
    Quad& quadrature;
    ArcLengthFunction arcLength;
    NextImpactFunction* nextImpact;
    NextImpactDiff* nextImpactDiff;
    Hunting* hunt;
    colvec initCond;
    double tolQuad, tolBis, tolNew, totalArc;
    int numberBounces;
    //@}

    /*! @name Parts of STA composite function. */
    //@{
    cube prepareReflectionMatrices();
    //@}
};

#endif// DYN_PLANARBILLIARD_H
